function osl_test_script_recon(testdatadir,testoutputdir,S)

global OSLDIR;

testoutputdir_full=[testoutputdir '_script_recon_' regexprep(S.forward_meg, ' ', '') '_hmm' num2str(S.do_hmm)];
runcmd(['rm -rf ' testoutputdir_full]);
mkdir(testoutputdir_full);

testplotsdir=[testoutputdir_full '/plots'];
mkdir(testplotsdir);

hmm_do_glm_statewise=1;
forward_meg=S.forward_meg;
do_hmm=S.do_hmm;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_beamformer_oat

printprefix='beamformer_wideband_oat';
printindex=1;

datadir=[testdatadir '/faces_subject1_data']; % directory where the data is

% Currently there is only 1 subject.
clear spm_files spm_files_epoched structural_files ;

% set up a list of mris
structural_files{1}=[datadir '/structurals/struct1.nii'];

% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[datadir '/spm8_meg1.mat'];
spm_files_epoched{1}=[datadir '/espm8_meg1.mat'];


%%%%%%
% DO REGISTRATION AND RUN FORWARD MODEL BASED ON STRUCTURAL SCANS


for i=1:length(spm_files),

    D=spm_eeg_load(spm_files{i});
    D = save_raw_tra_to_D(D);

    %D=spm_eeg_load(spm_files_epoched{i});
    %fidnew=D.fiducials;
    %fidnew.fid.label

    S2=[];

    S2.fid_label.nasion='Nasion';        
    S2.fid_label.lpa='LPA';
    S2.fid_label.rpa='RPA';

    S2.D = spm_files{i};    % requires .mat extension    
    %S2.D = spm_files_epoched{i};    % requires .mat extension    
    
    S2.mri=structural_files{i}; % set S2.sMRI=''; if there is no structural available        
    S2.useheadshape=1;

    %S2.sMRI=S2.mri;

    S2.forward_meg=forward_meg;

    S2.neuromag_planar_baseline_correction='vector_view';

    D=osl_forward_model(S2);
    
    
end;

ca;
%spm_eeg_inv_checkmeshes(D);
spm_eeg_inv_checkdatareg(D);
%spm_eeg_inv_checkforward(D, 1);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%
% SETUP BEAMFORMER AND FIRST-LEVEL GLM OAT
% In this section we will do a wholebrain beamformer, followed by a trial-wise
% GLM that will correspond to a comparison of the ERFs for the different
% conditions.

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.D_epoched=spm_files_epoched;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
oat.source_recon.freq_range=[1 48]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.4];
oat.source_recon.method=S.recon_method;
oat.source_recon.gridstep=12; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat.source_recon.forward_meg=forward_meg;
oat.source_recon.work_in_pca_subspace=S.work_in_pca_subspace;
oat.source_recon.hmm_epoch_after_prep=0;

if(do_hmm)
    oat.source_recon.hmm_num_states=-1;
    oat.source_recon.hmm_num_starts=2;
    oat.source_recon.hmm_pca_dim=30;
end;

design_matrix_summary={};design_matrix_summary{1}=[1 0 0 0];design_matrix_summary{2}=[0 1 0 0];design_matrix_summary{3}=[0 0 1 0];design_matrix_summary{4}=[0 0 0 1];
oat.first_level.design_matrix_summary=design_matrix_summary;

% Alternatively you can specify trial-wise subject-specific design matrices
% as ASCII files
if(0),
    % Need to create design matrix without rejects. D.pickconditions auto
    % excludes rejected trials:
    De=spm_eeg_load(spm_files_epoched{1});
    Ntrials=De.ntrials;
    X=zeros(Ntrials,length(oat.source_recon.conditions));
    tot=0;
    for ii=1:length(oat.source_recon.conditions),
        num=length(De.pickconditions(oat.source_recon.conditions{ii}));
        X(tot+1:tot+num,ii)=1;
        tot=tot+num;
    end;
    X=X(1:tot,:);
    
    % save as ASCII file:
    save([testoutputdir '/subject1_desmat.txt'],'X','-ascii');
    design_matrix_summary={};
    design_matrix_summary{1}=[testoutputdir '/subject1_desmat.txt'];
    
    try,rmfield(oat.first_level,'trial_rejects');catch, end;
    oat.first_level.design_matrix_summary=design_matrix_summary;
end;

% contrasts to be calculated:
oat.first_level.contrast={};
oat.first_level.contrast{1}=[3 0 0 0]'; % motorbikes
oat.first_level.contrast{2}=[0 1 1 1]'; % faces
oat.first_level.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
oat.first_level.contrast{4}=[0 0 -1 1]'; 
oat.first_level.contrast{5}=[0 -1 0 1]'; 
oat.first_level.contrast_name{1}='motorbikes';
oat.first_level.contrast_name{2}='faces';
oat.first_level.contrast_name{3}='faces-motorbikes';
oat.first_level.contrast_name{4}='fear-happy';
oat.first_level.contrast_name{5}='fear-neutral';
oat.first_level.cope_type='acope';
oat.first_level.hmm_do_glm_statewise=hmm_do_glm_statewise;

oat = osl_check_oat(oat);

oat.to_do=[1 1 0 0];
oat.first_level.name=['wholebrain_first_level'];
oat = osl_run_oat(oat);

if(do_hmm)
    figs=osl_oat_plot_hmm_states(oat);    
    for ff=1:length(figs),
        print(figs{ff}, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;
    end;
end;

%%%%%%
% OUTPUT SUBJECT'S NIFTII FILES

S2=[];
S2.time.reduce_time=1;
S2.resample_method='fft_ds'
S2.time_ds=2;
S2.time_range=[-0.1 0.3];
S2.data=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
[stats] = osl_reduce_data_to_visualize(S2);

S2=[];
S2.oat=oat;
S2.stats=stats;
S2.first_level_contrasts=[1,3]; % list of first level contrasts to output
[statsdir,times,count]=osl_save_nii_stats(S2);

%%%

con=1;
resamp_gridstep=2;
zslice=30;
vol=51;

figure;
map=ra([statsdir '/tstat' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol);
x1=squash(map,abs(map));
percfrom=96;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: tstat' num2str(con) ' vol' num2str(vol)];
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
map=ra([statsdir '/cope' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol);
x1=squash(map,abs(map));
percfrom=96;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: cope' num2str(con) ' vol' num2str(vol)];
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%

con=1;
zslice=30;
vol=80;

figure;
map=ra([statsdir '/tstat' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol);
x1=squash(map,abs(map));
percfrom=96;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: tstat' num2str(con) ' vol' num2str(vol)];
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
map=ra([statsdir '/cope' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol);
x1=squash(map,abs(map));
percfrom=96;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: cope' num2str(con) ' vol' num2str(vol)];
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%

con=3;
zslice=30;
vol=61;

figure;
map=ra([statsdir '/tstat' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol);
x1=squash(map,abs(map));
percfrom=96;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: tstat' num2str(con) ' vol' num2str(vol)];
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
map=ra([statsdir '/cope' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol); 
x1=squash(map,abs(map));
percfrom=96;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: cope d' num2str(con) ' vol' num2str(vol)];
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%

con=3;
zslice=30;
vol=1;

figure;
map=ra([statsdir '/pseudo_zstat_var_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
map=map(:,:,:,vol); 
x1=squash(map,abs(map));
percfrom=5;percto=90;
low=percentile((x1),percfrom);high=percentile((x1),percto);
overlay_act(flipud(squeeze(map(:,:,zslice))'), flipud(squeeze(bgmap(:,:,zslice))'),'red2yellow',0,[low high],[3000 8000]);
tit=['ERF: pseudo_zstat_var, vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,20:2:70));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%
% Investigating locations/regions of interest using a whole-brain OAT
% In this section we will interrogate the wholebrain OAT (run above) using 
% MNI coordinates.

oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat.first_level.name='wholebrain_first_level';
oat=osl_load_oat(oat.source_recon.dirname,oat.first_level.name);

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
% Set the mni_coord to a location of interest
mni_coord=[38,-64,-14]; 

% Find the nearest index of the beamformed voxels to the specified mni_coord 
dists=(sqrt(sum((stats.mni_coord-repmat(mni_coord,length(stats.mni_coord),1)).^2,2)));
[dist,vox_coord]=min(dists); % vox_coord is the voxel index of the beamformed voxels

% Load in the stats structure outputted by OAT first level trial-wise GLM analysis.
%stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
  
% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our two contrasts of interest:
figure;subplot(1,2,1);
con=1; errorbar(stats.times, squeeze((stats.cope(vox_coord,:,con))),squeeze((stats.stdcope(vox_coord,:,con))),'r'); hold on;
con=2; errorbar(stats.times, squeeze((stats.cope(vox_coord,:,con))),squeeze((stats.stdcope(vox_coord,:,con))),'g');
con=3; errorbar(stats.times, squeeze((stats.cope(vox_coord,:,con))),squeeze((stats.stdcope(vox_coord,:,con))),'b');
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
plot4paper('time (s)',[oat.first_level.cope_type]);
title([oat.first_level.cope_type]);

% Plot the 1-tailed t-stats over time for each contrast:
subplot(1,2,2);
con=1; plot(stats.times, squeeze((stats.cope(vox_coord,:,con))./stats.stdcope(vox_coord,:,con)),'r','LineWidth',2); hold on;
con=2; plot(stats.times, squeeze((stats.cope(vox_coord,:,con))./stats.stdcope(vox_coord,:,con)),'g','LineWidth',2); 
con=3; plot(stats.times, squeeze((stats.cope(vox_coord,:,con))./stats.stdcope(vox_coord,:,con)),'b','LineWidth',2); 
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
plot4paper('time (s)','1-tailed t-stat');
title([oat.first_level.cope_type]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%
% Investigating locations/regions of interest using a ROI OAT
% This section will re-rerun just the first level ERF analysis using a ROI 
% in the temporal occiptial fusiform cortex

% Load the wholebrain OAT (which was run above), to make use of the settings 
% and source_recon results already in there
oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat.first_level.name='wholebrain_first_level';
oat=osl_load_oat(oat.source_recon.dirname,oat.first_level.name);

% Give the first level analysis a new name to avoid copying over whole
% brain analysis:
oat.first_level.name='roi_first_level'; 

% Provide the roi mask as a niftii file:
oat.first_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
 
oat = osl_check_oat(oat);

oat.first_level.cope_type='acope';

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% Spatially average
S2=[];
S2.oat=oat;
if(strcmp(oat.first_level.cope_type,'cope')),
    S2.do_resolve_sign_ambiguity=1; 
end;
S2.stats_fname=oat.first_level.results_fnames{1};
if(0),
    S2.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
else
    % note that you can use a single coordinate (or a list of coords) as
    % well):
    mni_coord=[38,-64,-14]; 
    S2.mask_fname=osl_mnicoords2mnimask(mni_coord,1,'mask');
end;    
[stats,times,mni_coords_used]=osl_output_roi_stats(S2);

% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our contrasts of interest:
figure; subplot(1,2,1);
con=1; errorbar(stats.times, squeeze((stats.cope(1,:,con))),squeeze((stats.stdcope(1,:,con))),'r'); hold on;
con=2; errorbar(stats.times, squeeze((stats.cope(1,:,con))),squeeze((stats.stdcope(1,:,con))),'g');
con=3; errorbar(stats.times, squeeze((stats.cope(1,:,con))),squeeze((stats.stdcope(1,:,con))),'b');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
plot4paper('time (s)',[oat.first_level.cope_type]);
title([oat.first_level.cope_type]);
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');

% Plot the 1-tailed t-stats over time for each contrast:
% Note that these t-stats will be a bit inflated due to the pooling over
% voxels not taking into account the spatial correlation. However, strict 
% statistical validity is not a major concern at the first level; we worry
% about this instead at the group level (and use permutations stats)
subplot(1,2,2);
con=1; plot(stats.times, squeeze((stats.cope(1,:,con))./stats.stdcope(1,:,con)),'r','LineWidth',2); hold on;
con=2; plot(stats.times, squeeze((stats.cope(1,:,con))./stats.stdcope(1,:,con)),'g','LineWidth',2); 
con=3; plot(stats.times, squeeze((stats.cope(1,:,con))./stats.stdcope(1,:,con)),'b','LineWidth',2); 
plot4paper('time (s)','1-tailed t-stat');
a=axis; axis([-0.2 0.4 a(3) a(4)]);
title([oat.first_level.cope_type]);
legend('Motorbikes','Faces','Faces-Motorbikes','Location','NorthWest');

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%
%% ROI Time-frequency analysis
% This section will re-rerun the first level analysis using a ROI in 
% the temporal occiptial fusiform cortex, and do a time-frequency
% trial-wise GLM first-level analysis

% Load the wholebrain OAT (which was run above), to make use of the settings 
% and source_recon results already in there

oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat.first_level.name='wholebrain_first_level';
oat=osl_load_oat(oat.source_recon.dirname,oat.first_level.name);

% res=osl_load_oat_results(oat,oat.source_recon.results_fname{1});

% Give the first level analysis a new name to avoid copying over previous
% first-level analyses:
oat.first_level.name='roi_tf_first_level'; 
oat.first_level.bc=[1 1 0];

oat.first_level.tf_freq_range=[1 40]; % frequency range in Hz
oat.first_level.time_range=[-0.2 0.4];
oat.first_level.tf_num_freqs=12;
oat.first_level.tf_method='hilbert';
oat.first_level.space_average=0;
oat.first_level.tf_hilbert_freq_res=8;
oat.first_level.hmm_artifact_states=[];
oat.first_level.contrast{2}=[0 1 1 1]';
oat.first_level.contrast{3}=[-3 1 1 1]';
oat.first_level.contrast{4}=[0 1 0 0]';
oat.first_level.contrast{5}=[0 0 1 0]';
oat.first_level.contrast{6}=[0 0 0 1]';
oat.first_level.contrast{7}=[-1 1 0 0]';
oat.first_level.contrast{8}=[-1 0 1 0]';
oat.first_level.contrast{9}=[-1 0 0 1]';
oat.first_level.cope_type='cope';

oat.first_level.bc=[ 1 1 0 1 1 1 0 0 0];

if(1),
    oat.first_level.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
else
    % note that you can use a single coordinate (or a list of coords) as
    % well):
    mni_coord=[30 -66 -24]; 
    oat.first_level.mask_fname=osl_mnicoords2mnimask(mni_coord,1,'mask');
end;    
oat.first_level.hmm_do_glm_statewise=hmm_do_glm_statewise;

oat = osl_check_oat(oat);

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);

%%%%%%
% Spatially average the results over the ROI and plot the results

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% Spatially average
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.mask_fname=[OSLDIR '/std_masks/Right_Temporal_Occipital_Fusiform_Cortex'];
   
[stats,times,mni_coords_used]=osl_output_roi_stats(S2);

% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our two contrasts of interest:
figure;
con=2; 
imagesc(1:length(stats.times), stats.frequencies, squeeze((stats.cope(1,:,con,:)))');axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['cope' num2str(con)]);
figure;
imagesc(1:length(stats.times), stats.frequencies, squeeze(stats.cope(1,:,con,:)./stats.stdcope(1,:,con,:))',[0 20]);axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['tstat' num2str(con)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%

% Plot the time course of a single freq bin:
figure;
freqbin=nearest(stats.frequencies,8);
con=1; plot(stats.times, squeeze((stats.cope(1,:,con,freqbin))./stats.stdcope(1,:,con,freqbin)),'r','LineWidth',2); hold on;
con=2; plot(stats.times, squeeze((stats.cope(1,:,con,freqbin))./stats.stdcope(1,:,con,freqbin)),'b','LineWidth',2); 
con=3; plot(stats.times, squeeze((stats.cope(1,:,con,freqbin))./stats.stdcope(1,:,con,freqbin)),'g','LineWidth',2); 
legend('Motorbikes','Faces','Faces-Motorbikes');
xlabel('time (s)'); ylabel('1-tailed t-stat'); title([num2str(stats.frequencies(freqbin)) ' Hz']);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

if(do_hmm),
    for kk=1:5,
        figure;
        freqbin=nearest(stats.frequencies,8);
        con=1; plot(stats.times, squeeze(mean(stats.cope_by_state(1,:,con,freqbin,kk)./stats.stdcope_by_state(1,:,con,freqbin,kk),1)),'r','LineWidth',2); hold on;
        con=2; plot(stats.times, squeeze(mean(stats.cope_by_state(1,:,con,freqbin,kk)./stats.stdcope_by_state(1,:,con,freqbin,kk),1)),'b','LineWidth',2); 
        con=3; plot(stats.times, squeeze(mean(stats.cope_by_state(1,:,con,freqbin,kk)./stats.stdcope_by_state(1,:,con,freqbin,kk),1)),'g','LineWidth',2); 
        %con=1; plot(squeeze((stats.cope(1,:,con,freqbin))),'r','LineWidth',2); hold on;
        %con=2; plot(squeeze((stats.cope(1,:,con,freqbin))),'b','LineWidth',2); 
        %con=3; plot(squeeze((stats.cope(1,:,con,freqbin))),'g','LineWidth',2); 
        legend('Motorbikes','Faces','Faces-Motorbikes');
        xlabel('time (s)'); ylabel('1-tailed t-stat'); title([num2str(stats.frequencies(freqbin)) ' Hz']);
    end;
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% osl_example_continuous_oat

printprefix='beamformer_continuous_beta_oat';
printindex=1;

datadir=[testdatadir '/ctf_fingertap_subject1_data']; % directory where the data is

% Set up the list of subjects and their structural scans for the analysi 
% Currently there is only 1 subject.
clear spm_files structural_files;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files={[datadir '/dsubject1.mat']};
structural_files = {[datadir '/subject1_struct.nii']};      

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Undefined'};
oat.source_recon.freq_range=[13 30]; % frequency range in Hz
oat.source_recon.time_range=[300,32*30];
%oat.source_recon.time_range=[300,12*30];
oat.source_recon.method=S.recon_method;
oat.source_recon.gridstep=10; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.modalities = {'MEG'};
oat.source_recon.fid_label.nasion='nas';        
oat.source_recon.fid_label.lpa='lpa';
oat.source_recon.fid_label.rpa='rpa';
oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat.source_recon.forward_meg=forward_meg;
if(do_hmm)
    oat.source_recon.hmm_num_states=13;
    oat.source_recon.hmm_num_starts=2;
    oat.source_recon.hmm_pca_dim=30;
end;
oat.source_recon.mask_fname=[testdatadir '/top_brain_from40_2mm'];
oat.first_level.name='wholebrain';

oat = osl_check_oat(oat);

oat.to_do=[1 0 0 0];
oat.source_recon.pca_dim=100;
oat.source_recon.force_pca_dim=1;
oat.source_recon.regpc=0;
oat.source_recon.work_in_pca_subspace=S.work_in_pca_subspace;

oat = osl_run_oat(oat);

if(do_hmm)
    figs=osl_oat_plot_hmm_states(oat);    
    for ff=1:length(figs),
        print(figs{ff}, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;
    end;

end;

%%%
% Establish regressor for continuous time GLM. 

D=spm_eeg_load(spm_files{1});

% This should be setup to correspond to the same time window
% as the full time window for D

x=zeros(length(D.time),5);

block_length=30; %s
block_order=[5 5 5 5 5 5 5 5 5 5 4 3 2 1 2 3 1 4 3 4 1 3 2 1 4 4 2 1 3 3 4 1 4 3 1 2 1 2 3 4 3 4 1 2 3 4 1 2];

% [Left, Right, Rest, Both, Rest_at_start]
% [  1     2      4     8     16 ]
% figure;plot(D.time,squeeze(D(1,:,:)))
% emacs ~/vols_data/From_Nottingham_with_Love/JRH_MotorCon_20100429_01_FORMARK.ds/MarkerFile.mrk 

tres=1/(D.fsample);
tim=1;
for tt=1:length(block_order),    
    x(tim:tim+block_length/tres-1,block_order(tt))=1;
    tim=tim+block_length/tres;
end;

figure;plot(D.time,x);

print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%
% run glm to do regression against known finger tapping regressors

oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
oat.first_level.cope_type='cope';

% GLM stuff:
oat.first_level.design_matrix=x';
oat.first_level.contrast{1}=[-1 0 1 0 0]'; % rest-left
oat.first_level.contrast{2}=[0 -1 1 0 0]'; % rest-right
oat.first_level.contrast{3}=[0  0 1 -1 0]'; % rest-both
oat.first_level.tf_hilbert_freq_res=diff(oat.first_level.tf_freq_range);
oat.first_level.hmm_do_glm_statewise=hmm_do_glm_statewise;
oat.first_level.tf_method='hilbert';
oat.first_level.do_glm_demean=1;

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

% output niis
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1:3]; % list of first level contrasts to output
S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
[statsdir,times]=osl_save_nii_stats(S2);    

con=3;
figure;
resamp_gridstep=2;
map=ra([statsdir '/tstat' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=95;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
yslice=52;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000]);
tit=['Continuous: tstat' num2str(con) ' vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,40:2:80));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
con=3;
resamp_gridstep=2;
map=ra([statsdir '/cope' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=70;percto=99;
low=percentile((x1),percfrom);high=percentile((x1),percto);
yslice=51;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000]);
tit=['Continuous: cope' num2str(con) ' vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,40:2:80));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
con=3;
resamp_gridstep=2;
map=ra([statsdir '/pseudo_zstat_var_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=5;percto=90;
low=percentile((x1),percfrom);high=percentile((x1),percto);
yslice=51;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000]);
tit=['Continuous: pseudo_zstat_var, vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,40:2:80));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

% contrast_num=3;runcmd(['fslview ' statsdir '/tstat' num2str(contrast_num) '_2mm']);

%%%%
% Look at ROI time courses
% calculate an ROI mask 

oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];

% use FSL maths to threshold to create mask
con=3;
map=ra([statsdir '/tstat' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=95;percto=99.9;
low=percentile((x1),percfrom);high=percentile((x1),percto);
thresh=(high+low)/2;
runcmd(['fslmaths ' statsdir '/tstat' num2str(con) '_2mm -thr ' num2str(thresh) ' ' statsdir '/tstat' num2str(con) '_2mm_mask']);

figure;
map=ra([statsdir '/tstat' num2str(con) '_2mm_mask']);
resamp_gridstep=2;bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
low=0;high=high;
yslice=51;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000]);
title(['ROI mask, min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%
% look at raw timecourse averaged over ROI

oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];

source_recon_results=osl_load_oat_results(oat,oat.source_recon.results_fnames{1});
results = osl_get_recon_timecourses( source_recon_results, [statsdir '/tstat' num2str(con) '_2mm_mask'] );
time_ind=intersect(find(D.time>=oat.source_recon.time_range(1)),find(D.time<=oat.source_recon.time_range(2)));

figure;plot(results.times,normalise(squeeze(mean(results.source_timecourses(1,:),1))));
% compare to design matrix:
ho;plot(D.time(time_ind),x(time_ind,1),'r');
plot(D.time(time_ind),x(time_ind,2),'g');
plot(D.time(time_ind),x(time_ind,4),'k');
legend('data','left','right','both');
plot4paper('time(secs)','beta power');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%
% look at Hilbert Envelope timecourse averaged over ROI

oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];

oat.first_level.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];
oat.first_level.doGLM=0; % does not fit GLM and will output timeseries that would have been input into GLM
oat.first_level.tf_downsample_factor=1; 
oat.first_level.time_moving_av_win_size=4;
oat.first_level.space_average=0;
oat.first_level.name='roi';

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
figure;plot(stats.glm_input_times,normalise(squeeze(mean(stats.glm_input_data,1))));
% compare to design matrix:
ho;plot(D.time(time_ind),x(time_ind,1),'r');
plot(D.time(time_ind),x(time_ind,2),'g');
plot(D.time(time_ind),x(time_ind,4),'k');
legend('data','left','right','both');
plot4paper('time(secs)','beta power');
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

%%%%%%%%%
% Look at ROI wideband time-frequency spectogram
% we first will need to redo beamformer with a wider band

if(~do_hmm)
    
    oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
    oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
    statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
    con=3;
    oat.source_recon.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];

    printprefix='beamformer_continuous_wideband_oat';
    printindex=1;

    % need new OAT dirname for new source recon
    oat.source_recon.dirname=[testoutputdir_full '/' printprefix];

    oat.source_recon.freq_range=[4 48];
    oat.to_do=[1 0 0 0];
    oat = osl_run_oat(oat);

    %%%%%
    % now do first level for multiple freqs

    oat.first_level.mask_fname=oat.source_recon.mask_fname;
    oat.first_level.tf_method='hilbert';
    %oat.first_level.tf_morlet_factor=7;
    oat.first_level.tf_num_freqs=20;
    oat.first_level.tf_freq_range=[4 40];
    oat.first_level.tf_hilbert_freq_res=5;
    oat.first_level.doGLM=0; % does not fit GLM and will output timeseries that would have been input into GLM
    oat.first_level.space_average=0;
    oat.first_level.hmm_do_glm_statewise=hmm_do_glm_statewise;

    oat.to_do=[0 1 0 0];
    oat = osl_run_oat(oat);
    stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

    figure;imagesc(stats.glm_input_times,stats.glm_input_frequencies,squeeze(mean(stats.glm_input_data,1))');colorbar;
    axis xy;

    % compare to design matrix:
    ho;
    x3=x(:,3);
    x3(x3==0)=nan;
    plot(D.time(time_ind),10*x3(time_ind),'k','LineWidth',5);ho;
    legend('rest');
    plot4paper('beta power','time(secs)');
    print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

end;

%%%
%% DO seed based correlation

printprefix='beamformer_continuous_beta_oat';
printindex=1;
oat.source_recon.dirname=[testoutputdir_full '/' printprefix];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain');


printprefix='beamformer_continuous_beta_vector_oat';
printindex=1;
oat.source_recon.dirname=[testoutputdir_full '/' printprefix];

oat.source_recon.type='vector';


oat.first_level.name=['seedbased'];
oat.first_level.space_average=0;
oat.first_level.time_moving_av_win_size=2;%sec
oat.first_level.cope_type='cope';
oat.first_level.doGLM=1;
oat.first_level.tf_num_freqs=1;
oat.first_level.do_glm_demean=1;

% remove unwanted fields
try,oat.first_level=rmfield(oat.first_level, 'design_matrix');catch,end;
try;oat.first_level=rmfield(oat.first_level, 'contrast');catch,end;
try;oat.first_level=rmfield(oat.first_level, 'mask_fname');catch,end;

% seed info
%oat.first_level.connectivity_seed_mni_coord=[30 -30 48];
% seed info
oat.first_level.connectivity_seed_mni_coord=[42 -24 48];

oat.first_level.connectivity_seed_regress_zerolag=1;
oat.first_level.hmm_do_glm_statewise=0;
oat.first_level.connectivity_fc_methods='he';
oat.first_level.connectivity_seed_search_direction=0;

oat.to_do=[1 1 0 0];

oat = osl_run_oat(oat);

% output niis
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1]; % list of first level contrasts to output
S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
S2.resamp_gridstep=2;
[statsdir,times]=osl_save_nii_stats(S2);    

figure;
con=1;
resamp_gridstep=2;
map=ra([statsdir '/tstat' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=70;percto=99;
low=percentile((x1),percfrom);high=percentile((x1),percto);
yslice=51;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000],'jet',[20 high]);

tit=['Seed: tstat' num2str(con) ' vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,40:2:80));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;
con=1;
resamp_gridstep=2;
map=ra([statsdir '/cope' num2str(con) '_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=70;percto=99;
low=percentile((x1),percfrom);high=percentile((x1),percto);
yslice=51;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000],'jet',[low high]);
tit=['Seed: cope' num2str(con) ' vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,40:2:80));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

if(0)
figure;
con=1;
resamp_gridstep=2;
map=ra([statsdir '/pseudo_zstat_var_' num2str(resamp_gridstep) 'mm']);
bgmap=ra([OSLDIR '/std_masks/MNI152_T1_' num2str(resamp_gridstep) 'mm_brain']);
x1=squash(map,abs(map));
percfrom=5;percto=90;
low=percentile((x1),percfrom);high=percentile((x1),percto);
yslice=51;vol=1;
overlay_act(flipud(squeeze(map(:,yslice,:,vol))'), flipud(squeeze(bgmap(:,yslice,:))'),'red2yellow',0,[low high],[3000 8000]);
tit=['Seed: pseudo_zstat_var, vol' num2str(vol)]
title([tit ', min=' num2str(low) ', max=' num2str(high)]);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;

figure;lightbox(map(:,:,40:2:80));colorbar;title(tit);
print(gcf, '-dpng', [testplotsdir '/' printprefix num2str(printindex)]);printindex=printindex+1;
end;
% contrast_num=1;runcmd(['fslview ' statsdir '/tstat' num2str(contrast_num) '_2mm &']);


disp('*************************************************');
disp('Finished osl_test_scrip_recon test');
disp('*************************************************');